function [Ytest d out C] = fgmdmpval(COVtest,COVtrain,Ytrain,varargin)
    
    if isempty(varargin)
        method_mean = 'riemann';
        method_dist = 'riemann';
    else
        method_mean = varargin{1};
        method_dist = varargin{2};
    end
    labels = unique(Ytrain);
    Nclass = length(labels);
    
    % geodesic filtering
    [W,Cg] = fgda(COVtrain,Ytrain,method_mean,{},'shcov',{});
    COVtrain = geodesic_filter(COVtrain,Cg,W(:,1:Nclass-1));
    COVtest = geodesic_filter(COVtest,Cg,W(:,1:Nclass-1));    
    
    [Ytest, d, C] = mdm(COVtest,COVtrain,Ytrain,method_mean,method_dist);
    
    % estimation of center
    for i=1:Nclass
        C{i} = mean_covariances(COVtrain(:,:,Ytrain==labels(i)),method_mean);
    end

    % calcul distance train
    NTraintrial = size(COVtrain,3);
    dtrain = zeros(NTraintrial,Nclass);
    for j=1:NTraintrial
        for i=1:Nclass
            dtrain(j,i) = distance(COVtrain(:,:,j),C{i},method_dist);
        end
    end
    dtrain = log(dtrain);
    
    for i=1:Nclass
        mu(i) = mean(dtrain(Ytrain==labels(i),i));
        sd(i) = std(dtrain(Ytrain==labels(i),i));
    end
    
    % classification
    NTesttrial = size(COVtest,3);
    
    %for each sample we have Nclass distances (to choose from)
    d = zeros(NTesttrial,Nclass);
    for j=1:NTesttrial
        for i=1:Nclass
            d(j,i) = distance(COVtest(:,:,j),C{i},method_dist);
        end
    end
    
    d = log(d);
    
%     for i=1:Nclass
%         out(:,i) = 1-normcdf(d(:,i),mu(i),sd(i));
%         d(out<0.20) = NaN;
%     end
    
    for j=1:NTesttrial
        for i=1:Nclass
            
            %for each sample we know the distance to each class
            %now we check which distance fits well the class distribution
            p = normcdf(d(j,i),mu(i),sd(i));
            if ( (1 - p) < 0.10) % so we remove for sample j all the distances to classes which have low probability (the distances do not fit well to the class)
                 d(j,i) = NaN;
            end;
        end
    end
    
    [~,ix] = min(d,[],2);%NaNs are not considered
    Ytest = labels(ix);
    
    % nan class
    nc = sum(isnan(d'))==size(d,2);
    Ytest(nc) = NaN;
    
